%Create 
FileNamesOuter=cell(1,6);
FileNamesOuter{1,1} = {'output\FirstWorkspace6AH20','output\SecondWorkspace6AM20','output\ThirdWorkspace6AU20','output\FourthWorkspace6AZ20','output\FifthWorkspace6AH21'};
FileNamesOuter{1,2} = {'output\FirstWorkspace6BH20','output\SecondWorkspace6BM20','output\ThirdWorkspace6BU20','output\FourthWorkspace6BZ20','output\FifthWorkspace6BH21'};
FileNamesOuter{1,3} = {'output\FirstWorkspace6CH20','output\SecondWorkspace6CM20','output\ThirdWorkspace6CU20','output\FourthWorkspace6CZ20','output\FifthWorkspace6CH21'};
FileNamesOuter{1,4} = {'output\FirstWorkspace6EH20','output\SecondWorkspace6EM20','output\ThirdWorkspace6EU20','output\FourthWorkspace6EZ20','output\FifthWorkspace6EH21'};
FileNamesOuter{1,5} = {'output\FirstWorkspace6JH20','output\SecondWorkspace6JM20','output\ThirdWorkspace6JU20','output\FourthWorkspace6JZ20','output\FifthWorkspace6JH21'};    
FileNamesOuter{1,6} = {'output\FirstWorkspace6SH20','output\SecondWorkspace6SM20','output\ThirdWorkspace6SU20','output\FourthWorkspace6SZ20','output\FifthWorkspace6SH21'};
%Create string with root letters for different markets, A=AUD/USD,
%B=GBP/USD, C=CAD/USD, E=EUR/USD, J=JPY/USD,S=CHF/USD
Strcatroots = {'A','B','C','E','J','S'};

%Change directories
%functionsFolder=cd('..');
%parentDirectory=cd('output\');


nVec = [1 5 10 50 200 1000 5000];
nVecCounter = 0;

%Loop over different durations of regression
while(nVecCounter<size(nVec,2))
nVecCounter = nVecCounter+1;  
outerloop = 0;
%outerloop cycles through different roots (contract markets)
%can change outerloop to 6 from 5 to run CHF/USD test
while(outerloop<5)
    tic
    outerloop = outerloop + 1;
    FileNames = FileNamesOuter{1,outerloop};
    n = nVec(nVecCounter);

    %Initialize values for each storage variable
    
    %Beta Coefficients for different models and standard errors
    BetasMech = cell(11,1);
    BetasMech05d = cell(11,1);
    BetasNMech = cell(11,1);
    StandardErrorsMech = cell(11,1);
    StandardErrorsMech05d = cell(11,1);
    StandardErrorsNMech = cell(11,1);
    
    %Define general nonlinear model to estimate
    nlModel = @(Beta,X)(Beta(1)+Beta(2)*((abs(X).^1).*sign(X))); 
    %Define square root price impact model
    nlModel2 = @(Beta,X)(Beta(1)+Beta(2)*((abs(X).^0.5).*sign(X))); 

    
    LambdasNMechL = cell(11,1);
    BetasNMechL = cell(11,1);
    SandwichesNMechL = cell(11,1);
    DaySizesNMechL = cell(11,1);
    R2sNMechL = cell(11,1);

    LambdasMechL = cell(11,1);
    BetasMechL = cell(11,1);
    SandwichesMechL = cell(11,1);
    DaySizesMechL = cell(11,1);
    R2sMechL = cell(11,1);
    VolumesMechL = cell(11,1);
    AdditionalValues = cell(11,1);

    DatesForLater = cell(11,1);

    %Four vectors to build out: price with mechanical depletion, price with
    %only non-mechanical, signed order size, date/time (Y1Mech Y1NMech X1 t2)
    Y1MechAll = [];
    Y1NMechAll = [];
    X1All = [];
    t2All = [];
    retsAll = [];

    filecounter = 0;
    %Cycle through each contract month per market: Mar '20, June '20, Sep
    %'20, Dec '20, Mar '21
    while(filecounter<5)
        filecounter = filecounter + 1;
        clear ESU1820180801;
        
        %Load table of data for contract month
        load(FileNames{1,filecounter});
        
        %Convert timezome from UTC to NYC time
        a2 = datetime(ESU1820180801{:,1},'TimeZone','UTC');
        a2.TimeZone = 'America/New_York';
        t2 = datenum(a2);
        
        %Clean up NaNs if any
        dataIsNaNFormat = double(ESU1820180801{:,2:7});
        dataIsNaNFormat(dataIsNaNFormat == 0)= NaN;
        
        %Run routine that analyzes dataset to create a list of aggressive
        %orders signed by direction, with size and order book dynamics
        %corresponding to order time
        [AggressiveOrders2, ~,Y2] = MicroStructureES([t2 dataIsNaNFormat]);
        
        %Add a column that consists of the "true price"
        AggressiveOrders2(:,13) = (AggressiveOrders2(:,8).*AggressiveOrders2(:,9) + AggressiveOrders2(:,10).*AggressiveOrders2(:,7))./(AggressiveOrders2(:,8) + AggressiveOrders2(:,10));
        
        %Filter out trades outside of liquid trading hours
        t2H = hour(AggressiveOrders2(:,1));
        t2M = minute(AggressiveOrders2(:,1));
        
        selector = ((t2H>=18) +(t2H<15)); %FX & US Rates
        %selector = ((t2H>=3) +(t2H<12)); %Europe Rates
        AggressiveOrders3 = AggressiveOrders2(selector>0.01,:);
        clear AggressiveOrders2;
        
        %Create X and Y vectors using careful elimination of intraday jumps
        t2H = hour(AggressiveOrders3(:,1));
        diffjumps = ((abs(t2H(2:end)-t2H(1:end-1))>2)+ (abs(t2H(2:end)-t2H(1:end-1))<22))>1.01;
        diffjumps = [diffjumps; 0]; %#ok<AGROW> %next order is following day
        
        %Determine number of days in the dataset
        selectorDays = diffjumps.*(1:size(diffjumps,1))';
        selectorDays = [0;selectorDays(selectorDays>0,:);size(AggressiveOrders3,1)];
        
        %Pre-allocate based on number of days for daily estimates
        betasMech = zeros(size(selectorDays,1)-1,2);
        betasMech05d = zeros(size(selectorDays,1)-1,2);
        betasNMech = zeros(size(selectorDays,1)-1,2);
        standardErrorsMech = zeros(size(selectorDays,1)-1,2);
        standardErrorsMech05d = zeros(size(selectorDays,1)-1,2);
        standardErrorsNMech = zeros(size(selectorDays,1)-1,2);
        lambdasNMech = zeros(size(selectorDays,1)-1,1);
        betasNMechL = zeros(size(selectorDays,1)-1,2);
        sandwichesNMech = zeros(size(selectorDays,1)-1,4);
        daysizesNMech = zeros(size(selectorDays,1)-1,1);
        r2sNMech = zeros(size(selectorDays,1)-1,1);

        lambdasMechL = zeros(size(selectorDays,1)-1,1);
        betasMechL = zeros(size(selectorDays,1)-1,2);
        sandwichesMechL = zeros(size(selectorDays,1)-1,4);
        daysizesMechL = zeros(size(selectorDays,1)-1,1);
        r2sMechL = zeros(size(selectorDays,1)-1,1);
        volumesMechL = zeros(size(selectorDays,1)-1,1);
        additionalValues = zeros(size(selectorDays,1)-1,6);
        
        datesForLater = zeros(size(selectorDays,1)-1,1);
        
        counter = 0;
        while(counter<size(selectorDays,1)-1)
            %Loop through days and run analysis
            counter = counter + 1;
            startval = selectorDays(counter,1)+1;
            endval = selectorDays(counter+1,1);
            datesForLater(counter,1) = AggressiveOrders3(startval+1,1);
            %n=1; %n controls the number of transactions ahead for which we
            %are calculating price impact
            %X1 is signed order size
            X1 = AggressiveOrders3(startval:endval-n,3).*AggressiveOrders3(startval:endval-n,4);
            %Y1Mech is the difference in true price from before one trade
            %to the subsequent trade
            Y1Mech = AggressiveOrders3(startval+n:endval,6)-AggressiveOrders3(startval:endval-n,6);
            %Y1NMech is the non-mechanical change in true price (price
            %immediately before execution minus the price immediately after execution of the previous aggressive order
            Y1NMech = AggressiveOrders3(startval+n:endval,6)-AggressiveOrders3(startval:endval-n,13);
            %First with mechanical depletion
            try
                [beta,R,J,~,~,~] = nlinfit(X1,Y1Mech,nlModel,[0;0.01]);
                CovMtrxMech = NLLSVarMtrxCalcs(J,R);
                CovMtrxMech = CovMtrxMech/size(R,1);
                betasMech(counter,:) = beta';
                standardErrorsMech(counter,:) = (diag(CovMtrxMech).^0.5)';
            catch
                CovMtrxMech = NLLSVarMtrxCalcs(J,R);
                CovMtrxMech = CovMtrxMech/size(R,1);
                betasMech(counter,:) = beta';
                standardErrorsMech(counter,:) = (diag(CovMtrxMech).^0.5)';
            end
            %Now square root estimate
            try
                [betaMech05d,R105d,J05d,~,MSEMech05d,~] = nlinfit(X1,Y1Mech,nlModel2,[0;1]);
                CovMtrxMech05d = NLLSVarMtrxCalcs(J05d,R105d);
                CovMtrxMech05d = CovMtrxMech05d/size(R105d,1);
                betasMech05d(counter,:) = betaMech05d';
                standardErrorsMech05d(counter,:) = (diag(CovMtrxMech05d).^0.5)';
                betaMech05d'
            catch
                CovMtrxMech05d = NLLSVarMtrxCalcs(J05d,R105d);
                CovMtrxMech05d = CovMtrxMech05d/size(R105d,1);
                betasMech05d(counter,:) = betaMech05d';
                standardErrorsMech05d(counter,:) = (diag(CovMtrxMech05d).^0.5)';
                betaMech05d'
            end

           
        end
        
        %Store values
        Y1MechAll = [Y1MechAll;AggressiveOrders3(1+n:end,6)-AggressiveOrders3(1:end-n,6)]; %#ok<AGROW>
        Y1NMechAll = [Y1NMechAll;AggressiveOrders3(1+n:end,6)-AggressiveOrders3(1:end-n,13)];%#ok<AGROW>
        X1All = [X1All; AggressiveOrders3(1:end-n,3).*AggressiveOrders3(1:end-n,4)];%#ok<AGROW>
        t2All = [t2All; AggressiveOrders3(1:end-n,1)]; %#ok<AGROW>
        retsAll = [retsAll;(AggressiveOrders3(1+n:end,6)-AggressiveOrders3(1:end-n,6))./AggressiveOrders3(1:end-n,6)]; %#ok<AGROW>
        BetasMech{filecounter,1} = betasMech;
        BetasMech05d{filecounter,1} = betasMech05d;
        StandardErrorsMech{filecounter,1} = standardErrorsMech;
        StandardErrorsMech05d{filecounter,1} = standardErrorsMech05d;
        BetasNMech{filecounter,1} = betasNMech;
        StandardErrorsNMech{filecounter,1} = standardErrorsNMech;

        LambdasNMechL{filecounter,1} = lambdasNMech;
        BetasNMechL{filecounter,1} = betasNMechL;
        SandwichesNMechL{filecounter,1} = sandwichesNMech;
        DaySizesNMechL{filecounter,1} = daysizesNMech;
        R2sNMechL{filecounter,1} = r2sNMech;


        LambdasMechL{filecounter,1} = lambdasMechL;
        BetasMechL{filecounter,1} = betasMechL;
        SandwichesMechL{filecounter,1} = sandwichesMechL;
        DaySizesMechL{filecounter,1} = daysizesMechL;
        R2sMechL{filecounter,1} = r2sMechL;
        VolumesMechL{filecounter,1} = volumesMechL;
        AdditionalValues{filecounter,1} = additionalValues;
        
        DatesForLater{filecounter,1} = datesForLater;
    end

    %Having finished calculating values for all 5 contract months,
    %Cycle through and store values in vectors with daily estimates
    betasExtrasMech = zeros(300,2);
    stErExtrasMech = zeros(300,2);
    counter = 0;
    pointer = 1;
    while(counter<5)
    counter = counter +1;
    lengthBetas = size(BetasMech{counter,1},1);
    betasExtrasMech(pointer:pointer+lengthBetas-1,:) = BetasMech{counter,1};
    stErExtrasMech(pointer:pointer+lengthBetas-1,:) = StandardErrorsMech{counter,1};
    pointer = pointer + lengthBetas;
    end
    betasExtrasMech = betasExtrasMech(1:pointer-1,:);
    stErExtrasMech = stErExtrasMech(1:pointer-1,:);
    
    betasExtrasMech05d = zeros(300,2);
    stErExtrasMech05d = zeros(300,2);
    counter = 0;
    pointer = 1;
    while(counter<5)
    counter = counter +1;
    lengthBetas = size(BetasMech05d{counter,1},1);
    betasExtrasMech05d(pointer:pointer+lengthBetas-1,:) = BetasMech05d{counter,1};
    stErExtrasMech05d(pointer:pointer+lengthBetas-1,:) = StandardErrorsMech05d{counter,1};
    pointer = pointer + lengthBetas;
    end
    betasExtrasMech05d = betasExtrasMech05d(1:pointer-1,:);
    stErExtrasMech05d = stErExtrasMech05d(1:pointer-1,:);


    betasExtrasNMech = zeros(300,2);
    stErExtrasNMech = zeros(300,2);
    counter = 0;
    pointer = 1;
    while(counter<5)
    counter = counter +1;
    lengthBetas = size(BetasNMech{counter,1},1);
    betasExtrasNMech(pointer:pointer+lengthBetas-1,:) = BetasNMech{counter,1};
    stErExtrasNMech(pointer:pointer+lengthBetas-1,:) = StandardErrorsNMech{counter,1};
    pointer = pointer + lengthBetas;
    end
    betasExtrasNMech = betasExtrasNMech(1:pointer-1,:);
    stErExtrasNMech = stErExtrasNMech(1:pointer-1,:);


    lambdaExtrasNoMech = zeros(300,1);
    betaExtrasNoMech = zeros(300,2);
    sandwichesExtrasNoMech = zeros(300,4);
    daySizesNoMech = zeros(300,1);
    r2sNoMech = zeros(300,1);

    counter = 0;
    pointer = 1;
    while(counter<5)
    counter = counter +1;
    lengthLambda = size(LambdasNMechL{counter,1},1);
    lambdaExtrasNoMech(pointer:pointer+lengthLambda-1,1) = LambdasNMechL{counter,1};
    betaExtrasNoMech(pointer:pointer+lengthLambda-1,:) = BetasNMechL{counter,1};
    sandwichesExtrasNoMech(pointer:pointer+lengthLambda-1,:) = SandwichesNMechL{counter,1};
    daySizesNoMech(pointer:pointer+lengthLambda-1,:) = DaySizesNMechL{counter,1};
    r2sNoMech(pointer:pointer+lengthLambda-1,:) = R2sNMechL{counter,1};
    pointer = pointer + lengthLambda;
    end
    lambdaExtrasNoMech = lambdaExtrasNoMech(1:pointer-1,:);
    betaExtrasNoMech = betaExtrasNoMech(1:pointer-1,:);
    sandwichesExtrasNoMech = sandwichesExtrasNoMech(1:pointer-1,:);
    daySizesNoMech = daySizesNoMech(1:pointer-1,:);
    r2sNoMech = r2sNoMech(1:pointer-1,:);

    lambdaExtrasMech = zeros(300,1);
    betaExtrasMech = zeros(300,2);
    sandwichesExtrasMech = zeros(300,4);
    daySizesMech = zeros(300,1);
    r2sMech = zeros(300,1);
    volumesMech = zeros(300,1);

    counter = 0;
    pointer = 1;
    while(counter<5)
        counter = counter +1;
        lengthLambda = size(LambdasMechL{counter,1},1);
        lambdaExtrasMech(pointer:pointer+lengthLambda-1,1) = LambdasMechL{counter,1};
        betaExtrasMech(pointer:pointer+lengthLambda-1,:) = BetasMechL{counter,1};
        sandwichesExtrasMech(pointer:pointer+lengthLambda-1,:) = SandwichesMechL{counter,1};
        daySizesMech(pointer:pointer+lengthLambda-1,:) = DaySizesMechL{counter,1};
        r2sMech(pointer:pointer+lengthLambda-1,:) = R2sMechL{counter,1};
        volumesMech(pointer:pointer+lengthLambda-1,:) = VolumesMechL{counter,1};
        pointer = pointer + lengthLambda;
    end

    lambdaExtrasMech = lambdaExtrasMech(1:pointer-1,:);
    betaExtrasMech = betaExtrasMech(1:pointer-1,:);
    sandwichesExtrasMech = sandwichesExtrasMech(1:pointer-1,:);
    daySizesMech = daySizesMech(1:pointer-1,:);
    r2sMech = r2sMech(1:pointer-1,:);
    volumesMech = volumesMech(1:pointer-1,:);

    clearvars a2 AggressiveOrders3 dataIsNaNFormat diffjumps ESU1820180801 RegMtrx2 selector t2H t2M ;
    %plot(betasExtrasMech(:,3))
    

    %Calculate composite non-linear model (full)
    [trimmedData,selector] = VectorTrimmer(hour(t2All),[t2All X1All Y1MechAll Y1NMechAll retsAll],n);
    [betaMechComposite,R1,J,~,MSEMech,~] = nlinfit(trimmedData(:,2),trimmedData(:,3),nlModel,[0;1]);
    CovMtrxMech = NLLSVarMtrxCalcs(J,R1);
    CovMtrxMech = CovMtrxMech/size(R1,1);
    betasMech = betaMechComposite';
    standardErrorsMechComposite = (diag(CovMtrxMech).^0.5)';
    betaMechComposite'

    %Square root estimation section
    [~,Y,BetaHat,~,ehat,Sandwich] = OLSGenericProcessor(trimmedData(:,2),trimmedData(:,3));
    zMechComposite = (BetaHat./diag(sqrt(Sandwich)))*sqrt(size(trimmedData(:,2),1));
    r2 = 1- sum(ehat.^2)/sum((Y-mean(Y)).^2);
    lambdasMechLComposite = BetaHat(2,1);
    betasMechLComposite = BetaHat';
    sandwichesMechLComposite = Sandwich(:)';
    r2sMechLComposite = r2;
    MSELinearCompositeMech = mean(ehat.^2);

    
    %save(str);
    toc


holidays = [11 28; 1 20; 2 17; 5 31; 7 3; 9 7];
selector = volumesMech>=100000;

%SquareRoot estimation again
nlModel2 = @(Beta,X)(Beta(1)+Beta(2)*((abs(X).^0.5).*sign(X))); %Define nonlinear model to estimate
[betaMechComposite05,R105,J05,~,MSEMech05,~] = nlinfit(trimmedData(:,2),trimmedData(:,3),nlModel2,[0;1]);
CovMtrxMech05 = NLLSVarMtrxCalcs(J05,R105);
CovMtrxMech05 = CovMtrxMech05/size(R105,1);
betasMech05 = betaMechComposite05';
standardErrorsMechComposite05 = (diag(CovMtrxMech05).^0.5)';
betaMechComposite05'

%Store values
datesExtras = zeros(300,1);
counter = 0;
pointer = 1;
while(counter<5)
counter = counter +1;
lengthdates = size(DatesForLater{counter,1},1);
datesExtras(pointer:pointer+lengthdates-1,:) = DatesForLater{counter,1};
pointer = pointer + lengthdates;
end
datesExtras = datesExtras(1:pointer-1,:);

%Filter out holidays
selectorNonHolidays = zeros(315,1)+1;
selectorNonHolidays([8 63 127 192 257],:) = 0;
selectorNonHolidays = (selectorNonHolidays>0);
formatOut = 'mm/dd/yy';
dateVec = datestr(datesExtras+1,formatOut);

str = strcat('output\','ComprehensiveStudy_6',Strcatroots{1,outerloop},'All',num2str(nVec(1,nVecCounter)),'OLSOnly');

%Save to file
save(str);

end

end

%Redirect back to main directory
%cd(parentDirectory);